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The sputtering yield, Y , from a cylindrical thermal spike is calculated using a two dimensional 
fluid dynamics model which includes the transport of energy, momentum and mass. The results show 
that the high pressure built-up within the spike causes the hot core to perform a rapid expansion 
both laterally and upwards. This expansion appears to play a significant role in the sputtering 
process. It is responsible for the ejection of mass from the surface and causes fast cooling of the 
cascade. The competition between these effects accounts for the nearly linear dependence of Y with 
the deposited energy per unit depth that was observed in recent Molecular Dynamics simulations. 
Based on this we describe the conditions for attaining a linear yield at high excitation densities and 
give a simple model for this yield. 

PACS numbers: 79.20.Rf, 47.40. Nm, 83.85.Pt 
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I. INTRODUCTION 



The ejection of atoms from the surface of a solid during 
ion irradiation ismwell documented both experimentally 
and theoreticallytl. This phenomenon, known as sput- 
tering, is due to the energy transferred to the atoms in 
the target by the incident ion. This produces a cascade 
which can cause some atoms to overcome the surface's 
attractive barrier and escape to vacuum. 

In previous theoretical work the mean number of 
ejected atoms per incoming ion, or sputtering yield Y , 
is related to the energy deposited the ion per unit thick- 
ness at the surface of the target Fd, as oc Fg. The 
value of the power n depends on the type of collision cas- 
cade produced by the ion, namely linear and non-linear 
cascades. For linear cascades, when the density of mov- 
ing atoms N„iov within the cascade is^ small compared 
to normal density Nq, one has n — laH, whereas in the 
non-linear case Nmov ~ theotetical work predicted 
that n must be greater than oneuQ. These results are so 
firmly established that the consensus among workers in 
the field is that n >_1 and non-linear cascades are to some 
extent synonymousO. Similar results have been found for 
sputteong in response to electronic energy deposited in 
a solidQ, but here we refer to work on collision cascade 
sputtering. 

Recent Molecular Dynamics (MD) studiesoO, however, 
cast doubt on this relationship. According to these pa- 
pers, purposely prepared non-linear cascades can give rise 
to sputtering yields which depend linearly oruJ^ d (see 
figures . Further evidence is found in Ref .113. After 
modifying the standard thermal spike theory (STST) to 
include the transport of mass, the sputtering yields so 
calculated appeared to be much closer to a linear func- 



tion of Fd than to the prejdicted by the STST. 

Although the results in Ref.t3 show the importance of 
having a target which can change its specific volume as a 
fluid, it is not a full fluid dynamics calculation. Since 
the target was assumed to be infinite, the sputtering 
yields had to be calculated in the same manner as in 
the STST. That is, an expression for the evaporation 
rate was used that was borrowed from the kinetic theory, 
and the sputtering yields were obtained by integrating it 
along a plane representing an otherwise non-existent sur- 
face. Further, the transport was only radial, but the MD 
calculations showed the importance of energy transport 
along the track. 

In order to circumvent such a difficulty, our previous 
calculations are extended to a target which, in addition 
to being compressible, has a solid-vacuum interface. To 
this end, the target density, velocity and internal energy 
are all assumed to vary with time in a manner which 
is described by the fluid dynamics equations. Conse- 
quently, sputtering emerges naturally, as that part of the 
target that succeeds in escaping from the condensed to 
the gaseous phase. 

The aim of this paper is to show the most relevant 
aspects of those calculations, from the underlying math- 
ematics to the results and implications of the proposed 
model. Although the present calculations can be applied 
to a variety of ion- induced thermal spike cases, we have 
purposely limited ourselves to the cases contained in pre- 
vious MD simulationsErEl. Therefore, the results in this 
paper are intended for cylindrical thermal spikes. In Sec- 
tion 1^ we introduce the fluid dynamics equations as well 
as the various expressions used along the present calcu- 
lations. Results and discussions are presented in Section 
III. Finally, the conclusions and suggestions for further 



studies are presented in Section IV 
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II. THEORY 

We assume that the target is a continuous medium 
with cyhndrical symmetry, and it is completely charac- 
terized by its atomic number density N, velocity v, and 
internal energy e (per atom) defined as, 



U 



(2.1) 



(2., 



The heat conduction coefficient is replaced by that in 
Ref.E 
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(2.9) 



where is the Boltzmann's coefficient, T the temper- 
ature and U is the potential energy per atom. It must 
be noted that by using the equation above the heat ca- 
pacity at constant volume, Cy , is assumed to be that of 
a dilute gas, i.e. 3fcB/2. This approximation however is 
fairly acceptaiJe for the purpose in this paper, since as 
shown in Ref.Ell, the quadratic dependence of Y with Fd 
does not appear to be connected to Cy ■ Moreover U is 
obtained from the expressionEj 
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(2.2) 



where M is the mass of the target atom, Cq is the speed of 
sound at T = O-ftT, and iVo is the normal atomic number 
density, /i and v are_two numerical constants which, as 
we explained in Ref.EJ, are not independent. Thus we set 
H — 2, then ly — + Mc^/Uq, Uq being the potential 
energy at normal density, i.e. Uq = pi7(7Vo). 

Using the same notation as in Ref.Ej, we write the fluid 
dynamics equations as follows 
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where the subscripts stand for the r and z copidinates, 
p is the pressure and cr^'j, is the viscosity tensorEj defined 
as. 



^\dxk dxi 



(2.6) 



where rj is the dynamic viscosity coefficient and Qcon and 
Qvis account for the heat produced by thermal conduc- 
tion and viscosity per unit volume and time, namely 



Qcon = V (ktVT) , 

where kt is the thermal conductivity and. 



(2.7) 



where 7ra = 1.151A2. This form was used in order to 
compare results and because there seems to be no reaspii 
for using a more "realistic" one since, as shown in Ref.EJ, 
kt and the quadratic dependence of the sputtering yield 
appear not to be connected. 

Making use of the fact that, for dilute gases, ry and kt 
are related through the equation 77 = k^M / {^ths), we 
thus introduce the dimensionless viscosity coefficient 



ri* = 3fcs77/(M«;T) 



(2.10) 



Similarly, the pressure p is assumed to be a function 
of both temperature pad density. Here, we follow the 
approximation in Ref.E^I and split p into two terms 



p = PT+Pc , 



(2.11) 



where the thermal pressure px can be obtained from the 
expression 



Pt = X NksT, 



(2.12) 



A being a numerical constant. The so called crystal 
pres sure pc can be obtained, from the potential energy 
Eq.(^.2[) using the equationll3 



PC 



, dU 

'on 



(2.13) 



For computational purposes, Eqs. ( |2.3[j2.5| ) are ap- 
plied to a finite system, which is defined by inequalities 
< r < Rmax and < z < z^ot (see figure |l|) . Further- 
more, the top wall, i.e. 2; = 0, is assumed to be made of 
a perfectly absorbent material, whereas the boundary at 
the bottom is perfectly closed as far as to the exchange of 
mass, momentum and energy is concerned. The lateral 
wall can be made either closed, like the bottom surface, 
or partially open. That is, closed to mass transport but 
open to energy and momentum exchange. Results in this 
paper were obtained using the latter option. Otherwise, 
in the case of a large deposited energy, one would need an 
exceedingly large target to minimize the effects of energy 
and momentum reflections. A more detailed-|description 
of this program will be published elsewherai^. 



2 



top boundary 









rarified fluid 


















































r 


ormal 


density fluid 






























































































































bottom wall 



Figure 1 ; Frame ot reference and grid. /Jakas et al./ 

FIG. 1. Sketch illustrating the frame of reference and grid 
utilized in the present calculations. At t=0 the "fluid" occu- 
pies the region defined by inequalities z^urf < z < Zf,ot and 
< i? < Rmax, and the hot spike is confined to a cylinder of 
radius Rcyi- 

At t = the target is within a range of z defined by in- 
equality (z > Zsurf)- For numerical reasons however, we 
assume that the region that would normally be a vacuum 
is filled with a low-density fluid, i.e. Nmin = x iVo. 
Exchange of energy, momentum and matter is forbidden 
in this fluid, as well as in any other piece of a fluid with 
density lower than A'niin- The possible net flux of matter 
is continuously checked along the fluid, and the restric- 
tions above are relaxed as soon as an element of the fluid 
would have its density increased above Nmin- 

To energize the spike, all atoms within a cylinder of 
radius Rcyi are given an energy E^xc above their initial 
energy, eo = —Uq + (3/2)fcBTo, where Tq is the back- 
ground temperature, often assumed to be 10 K. These 
are the ini 
simulation; 



conditions used in a number of the MD 
, again allowing direct comparison with the 
results here. The initial conditions for Eqs.(pT3-2.5) thus 
become. 



iV(0,r,z) 
e(0,r,z) 



No if 



Z ^ ^surf 

otherwise 



(2.14) 



Eexc + £0 if r < Rcyi and z > Zsurf 
U{N^,n) + (3/2)A:Bro if < z < 
£0 otherwise 



With the assumptions above, the deposited energy be- 
comes, 

Fd = TTRlyiNoEcxc ■ 



As is customary, in solving the fluid dynamics equa- 
tions the functions A^, v and e are defined over a discrete 
set of NR X NZ points, whose mesh size is determined 
by Ar and Az (See figure |l| and Table ||). A compromise 
has to be made about target size since a large target 
implies a fairly large system of coupled equations with 
fairly long running times. Whereas too small a target 
gives rise to boundary effects that would make calcula- 
tions meaningless. Similarly, in choosing Zsurf one has 
to take into account that during ejection not all the mat- 
ter that crosses the surface will be ejected. Therefore, 
the distance between the initial surface and the top wall 
should be large enough to not "collect" matter that, oth- 
erwise, would not be ejected. Finally, the piece of matter 
representing the target must be thick enough. The con- 
densed phase is assumed to be lOcr thick, which means 
that Zsurf -lOo- and Zbot «20(t. NR = 40 and NZ = 20 
were found to be adequate for all the cases studied in this 
paper. 

When in tegrating the fluid dynamics equations 
(p~3 2.4,2.5) from t —0 to tf, the total flux of matter 
passing through the top boundary is also calculated. In 
this way the sputtering yield is obtained as a function 
of time, Y{t). This is used to verify if i/ was long 
enough so that no matter remains within the system 
that may significantly contribute to the sputtering yield. 
We use the Y{tys for t < tf to extrapolate Y{t) to in- 
finity, i.e. Yoc ~ liiatf^oo Y{tf). Only runs for which 
Yoo — Y{tf) « O.lFoo are accepted. Normally, tf ranging 
from 10 up to 50 ps are required. 

Since calculations in this paper are meant to be com- 
pared with those in MD simulations, which often use 
Lennard-Jones (LJ) potentials, the various parameters 
characterizing our system correspond to those of Argon 
(see Table |). M^dflu and t/o=0.08eV have become stan- 
dard parametersEro although the LJ calculations fully 
scale with Uq and M. Therefore, the results apply to 
a broadet-set of materials as shown also using a Morse 
potentialtj. Consistent with this, for most cases we used 
Ar = Az — a, where a is the LJ distance. However, as 
several approximations were introduced, we cannot en- 
sure that the fluid in our calculations accurately describes 
the potentials used in the MD simulations. Similarly, we 
do not want to leave this section without mentioning that 
although the fluid re pres enting the target is assumed to 
be compressible, Eq.( |2.6| ) looks the same as that of an 
incompressible fluid because the Stokes' condition is as- 
sumed to hold. 

TABLE I. Value of the parameters used in the present cal- 
culations. 



Property 


Symbol 


Value 


Atomic mass 


M 


40.0 u.m.a. 


Atomic number density 


No 


0.026 at/A^ 


Speed of sound 


Co 


17 A/ps 


Binding energy 


Uo 


0.08 eV 


Lennard-Jones distance 


a 


3.405 A 
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III. RESULTS AND DISCUSSIONS 



We calculated the sputtering yield for different values 
of A, 77* and the speed of sound Co, and the results are 
depicted in figures || to ^. We observe that in all the 
cases the yield increases with increasing excitation en- 
ergy Eexc- Similarly, Eexc ~ Uq is an effective threshold 
for ejection for the initial radius used, since the yields 
rapidly decrease for Eexc comparable to or less than Uo- 
Whereas the MD requires varying potential types to ob- 
tain different material properties, here we do this by di- 
rectly varying the material properties. In this manner 
the relationship between different materials can be de- 
scribed. 

We observe that A has a great influence on the sput- 
tering yield. The larger the A the greater the yield. A=4 
appears to reproduce MD-simulations quite well, whereas 
A=2 and 1 underestimate the yields at small excitation 
energies. These results are, to some extent, easy to un- 
derstand: with all other parameters remaining the same, 
as A becomes larger the thermal pressure build up within 
the spike increases and more ejection is expected. 

It mu st be noted, however, that the total energy, c.f. 
Eq.( ^.l|) , does not depend on A. Increasing A only in- 
creases the thermal pressure and speeds up the con- 
version of thermal motion into directed kinetic energy. 
Therefore, thermal conductivity has less time to take en- 
ergy away from the spike and the ejection of matter in- 
creases. 





Figure 3: Sputtering yield for different viscosity coefficient (n ). /Jakas et al./ 

FIG. 3. Sputtering yield as a function of tlie excitation en- 
ergy and different values of viscosity coefficient ri* . 




Figure 4: Sputtering yield for different speed of sound. /Jakas et al./ 

FIG. 4. Sputtering yield as a function of the excitation en- 
ergy and different values of the speed of sound cq. 



Figure 2: Sputtering yield for different thermai pressure coefficient (1.). /Jakas et al./ 

FIG. 2. Sputtering yield as a function of the excitation en- 
ergy and different values of parameter A. 



The role played by the viscosity on the sputtering yield 
is illustrated in figure |. For the cases studied here, we 
observed that viscosity has a negative influence on the 
ejection process, as yields are seen to get smaller with 
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an increase of the viscosity coefficient. At small excita- 
tion energies the viscosity appears to play a major role. 
Furthermore, calculations using rj* = 0.1 produced a 
good agreement with MD-simulations while those with 
rj* = 0.2 and 0.4 resulted in significantly smaller yields. 
The fact that the best agreement with MD-simulations 
corresponds to calculations with rj* — 0.1 is not unex- 
pected since ry*-values of approximately ttift order have 
been calculated for a Lennard-Jones fluidli3. 



Finally, modifying the speed of sound does not pro- 
duce a significant change in the sputtering yield. Figure 
^ shows results for the speed of sound both above and be- 
low its normal value. The change in the sputtering yield 
is very small compared that produced by changing either 
the viscosity coefficient or the thermal pressure coeffi- 
cient A. We observe that, for high excitation energies, an 
increase in the speed of sound leads to a slightly greater 
yield, and that such a trend is reversed as Eexc/Uo be- 
comes smaller than 3. 
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Fig 5: Density' and mass-flin. /Jakas et al./ 

FIG. 5. These plots illustrate the density and mass-flux 
vectors at difi^erent times for a spike with dE /dX—4eY / K, 
A=4, ?7*=0.1, co=17 A/ps and Rcyi = 2a. The scale used for 
translating from relative density (N/Nq) into colors is shown 
up in the figure. Note that the scale is non-linear, as more 
colors are used at both small densities and around N/No — 1. 
Furthermore, due to interpolation in the plotting software, 
details of the order of the grid size, or smaller, might not 
be accurately copied. The horizontal line denotes the initial 
position of the surface. 



Fig 6: Temperature and niass-flnx. /Jakas el al./ 

FIG. 6. Temperature and mass-flux vectors at different 
times within the fluid for the spike described in caption ^. 
The scale used for color vs. temperature appears up in the 
figure. 
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Fig 7: Details of the ejection process. /Jakas et al./ 



FIG. 7. Close-ups of plots in figure q illustrating in more 
details the dynamics of the fiuid within the "core" of the spike 
and near the surface. 
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As we mentioned in the introduction, the most inter- 
esting result in this paper is our abihty to explore the ma- 
terial parameters that lead to the near linearity exhibited 
in the yield in our MD calculations even though the sput- 
tering is a non-linear process. By exploring the param- 
eter space we can better explain that phenomenon and 
assess its relevance. Our calculated yields in figures |[^ 
clearly show that a linear region is attained for Ee^c > C^o 
using a set of materials parameters. Therefore, non-linear 
sputtering does not necessarily imply non-linear yields. 
From these figures, it also appears that linearity is ap- 
proached at higher energy densities than those studied 
here for other materials parameters. Below we describe 
this phenomenon. 

To understand the change in the dependence of the 
yield with increasing excitation density for fixed Rcyi , we 
analyzed the time-evolution of the spike paying particular 
attention to those aspects of the energy and momentum 
transport that are related to the ejection of matter. To 
this end, in figures |^|| we have plotted the density, the 
mass-flrrx vector and temperature in the fluid at different 
times after the on-set of the spike. These cases corre- 
spond to a deposited energy of 4 eV/A, A = 4, r]*=0.1, 
co=17 A/ps and Rcyi=2a; and, in the three flgures, the 
initial surface is located at lOu, i.e. Zsurf = lOcr. 

One readily observes that the temperature within the 
spike drops below 500 K in approximately 1 ps, and 
that the fluid immediately surrounding the spike hardly 
reaches temperatures higher than, say, 100 K. This is 
in agreement with oup-MD results and our earlier fluid 
dynamics calculationsEJ. These studies already showed 
that, due to the quick, adiabatic expansion of the fluid, 
the temperature of the spike decreases much more rapidly 
than it would due to thermal conduction. In addition, for 
times greater than 1 ps the thermal energy appears to be 
converted into an clastic wave (seen as red in Fig. ||) 
that travels in the radial direction at approximately the 
speed of sound. The reader must be aware of the non- 
linear scale used in Fig. ^ where colors were purposely 
chosen so as to change rapidly around both A^o f^nd at 
low density. Due to this, even the rather small relaxation 
of the surface density appears as a yellow stripe, which 
extends to the right of the spike and gets thicker with in- 
creasing time. These figures show us that the whole pro- 
cess would be better described as an "explosion" rather 
than a smooth, therrrially diffusive release of energy as 
proposed in the STSTeI. 

Note that, in contrast to material further away from 
the surface, the fiuid that is near the surface and within 
the spike, appears to follow a spherical, rather than a 
cylindrical expansion. That is, if one interpolates the 
mass-flux vector and figures out the streamlines of the 
fluid, then, one can readily see that near the open bound- 
ary of the spike, they seem to radiate out from a point 
located on the spike axis somewhere below the surface. In 
order to understand this, one has to realize that the mo- 
mentum acquired by any particle within the fluid results 
from the fast, though small, displacements of the lateral 



and top boundaries which takes place at an earlier stage 
of the aforementioned explosion. 

The forces produced by such displacements propagate 
at the speed of sound which, within the hot spike, is faster 
than cocZl. Therefore, by the time all the fluid within the 
spike has been set into motion, i.e. t — Rcyi/c after the 
onset of the spike, a particle at (r, z) with 0< z < Rcyi 
and 0< r < Rcyi will have acquired a velocity that is 
proportional to the time it has been exposed to such 
forces, namely u,- oc r and Vz oc —{Rcyi — z). There- 
fore, as Vz/vr ~ —{Rcyi — z) j T this particle will appear 
as moving away from a point located exactly on the axis 
at Rcy/ below the surface. By the same token, any par- 
ticle at a depth greater than Rcyi within the spike, will 
remain unaware of the presence of the surface and its ve- 
locity will be directed along the radial direction (see fig- 
ure 0). With increasing time our description above will 
become less and less accurate. However, as the forces 
acting within the spike take the largest values during the 
earliest stage of the "explosion" , the velocities achieved 
by the fluid during that time essentially determine the 
subsequent dynamics of the spike. 

Another aspect of the velocity field which deserves at- 
tention is that around the rim, on the cold side of the 
spike. Contrary to what happens deep in the target, 
where the cold side is compressed and subsequently dis- 
placed along the radial direction, the rim is partially 
wiped out. This not only adds more matter to sputtering 
per se, but also clears the way for further ejection as it 
widens the sputtering radius, or the radius from within 
which particles are ejected. 
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FIG. 8. Schematics used to obtain the sputtering radius. 
From this simple picture one can readily calculate the 
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sputtering radius. To this end, we define the excess en- 
ergy as the total energy per particle relative to the bot- 
tom of the potential well, i.e. e — e + ^Mv^ + Uq. If 
one assumes that the elastic wave in the upper part of 
the spike propagates isentropically along the streamlines, 
one may write 



eA/d\ 



eB/d% 



(3.1) 



where and are the distance from the center of 
the spherical expansion to points A and B, respectively 
(see figure similarly, ba and bb are the correspond- 
ing excess energies. Therefore, as e >Uo is a necessary 
condition for ejection, ca = Eexc and Rcyi/dA^ Rs/dB, 
one can obtain the sputtering radius [Rs] aal3 




E_/U„ 



Rs ~ Rcyl {Eexc/Uo 



,1/2 



(3.2) 



Figure 9: Yields tor cjitferents Epil<e radii /Jakes et ai./ 



Accordingly, the sputtering yield can be calculated as 
the amount of mass contained within a cone of height 
Rcyl and base radius Rs, i.e. 



^cyl ' 



(3.3) 



In order to verify this simple expression, we calculate 
the sputtering yield for different spike radii. The results, 
that appear in figure ||, show that our fluid dynamics 
calculati ons c ompare fairly well with the MD yields, and 
that Eq.(3^) accounts reasonably well for the yields at 
high-excitation energies. Discrepancies between MD and 
fluid dynamics at low excitation energies and for small 
spike radii do appear. A detailed analysis of such devia- 
tions was not carried out. As previously mentioned, the 
various quantities entering our model do not accurately 
account for the Lennard-Jones fluid in the MD simula- 
tions. In addition, having assumed the solid target is 
a fluid, effects arising from the crystalline structure and 
the atomic nature of the target can not be described. In 
the MD simulations focused collision sequences play an 
important role at carrying energy and momentum away 
from the spike, particularly for sma ll sp ike radii. Finally, 
it is worth noticing that equation (3.3) predicts a linear 
dependence of the yield with the excitation or deposited 
energy. A result that was derived in Ref.Ej using a sim- 
ple, intuitive model rather than well supported, rigorous 
calculation. 



FIG. 9. Sputtering yield for different spike radii. MD cal- 
culations appear as symbols whereas Hydrodynamics results 
are plotted as thick lines. Thin, stra ight lines show the sput- 
tering yield obtained using Eq. (3.3). 



Although we have chosen not to address to the prob- 
lem of crater formation, it is worthwhile observing that 
late in our calculations craters do appear, and they are all 
surrounded by a rim several a high (the reader is referred 
to Ref.o for additional information about crater forma- 
tion) . The hole left by the spike is normally greater than 
the initial radius of the hot core. It is formed as the re- 
sult of the net displacement produced by the elastic wave 
along the radial direction. Near the edge of the hole, the 
radial momentum is less than it is in the material be- 
low. As a result, a kind of cantilever is formed which is 
pushed upwards by the fluid below. See the case of t=2 
ps in figure |^. It is important to mention, however, that 
for Eexc smaller than, say U, then no hole is formed. 



IV. CONCLUSIONS 



Sputtering at relatively high excitation densities is an 
old but unsolved theoretical problem in ion solids inter- 
actions. Analytic diffusive thermal spike models are com- 
monly used to interpret data at high excitation densities, 
although these models were never tested against more 
detailed calculations. In addition, there is a history of 
applying ideas from fluid dynamics to explain sputtering 
at high excitation density. These mcjuiels are called by_a. 
number of names (gas flowEiJ, shockllj, pressure pulseEHl 
etc.) and attempt to account for the fact that sputtering 
at high excitation density does not occur on an atom by 



7 



atom basis. These models also require a more detailed 
theoretical justification. 

Establishing a theoretical basis for sputtering models 
at high excitation density has been addressed recently 
by MD simulations using model materials and simplified 
descriptions of the initial conditions. It was shown that 
standard spike models break down at precisely those high 
excitation densities which they were intended to treat. In 
addition, a new sputtering regime was found. On increas- 
ing the energy density in the spike for fixed spike radius, 
the yield changed from a non-linear dependence on the 
excitation density to a linear dependence even though 
the ejection process clearly remained non-linear. This is 
contrary to the conventional wisdom and suggests satu- 
ration occurs in the sputtering. To examine this result 
we first showed that such a regime is never attained for 
any set of nmterial properties using the diffusive thermal 
spike modeltil. Since the standard spike model involves 
solving only the energy equation, we then numerically 
integrated the full se±-pf fluid equations for a ID model 
of a cylindrical spikellj. Differences with the MD result 
remained which we attributed to the lack of a surface. 
Here use a 2D fluid dynamics model with a surface to 
confirm that when the full set of equations is treated the 
MD result at high excitation density can be attained. 
Therefore, as pointed out earlier, the principal deficiency 
of the standard spike model is the assumption that the 
transport is diffusive. 

We have calculated for the first time the sputtering 
yield from a cylindrical thermal spike by directly inte- 
grating the full 2D fluid-dynamics equations. The trans- 
port of mass and momentum are seen to play a significant 
role in the ejection process. Since the conversion of the 
thermal energy into kinetic/potential energy within the 
spike occurs very early, the ejection process at high en- 
ergy densities is much more closely related to an "explo- 
sion" rathec-than to the thermal diffusion and evapora- 
tion modelsO typically used to describe sputtering at high 
energy densities. Comparisons with MD-simulations us- 
ing appropriate material parameters, show that our fluid 
dynamics description can account for the main features 
of the cylindrical thermal spike. These calculations also 
confirm the MD result that transport along the cylindri- 
cal axis is as important as radial transport and, therefore, 
a 2D model is required. We show the reported nearly 
linear yield comes about because of the competition be- 
tween pressurized ejection and the transport of energy 
away from the spike by the pressure pulse. 

Using the evolution of the streamlines seen in these 
calculations we obtain a simple expression for the yield 
at high excitation density for a reasonablc-set of mate- 
rial parameters. Bringa and co-workers (li3,Ej) had shown 
that in this regime the yield could be written in the form, 
Y « C[Rcyilir{[dEldx\eff{l/UQ)Y where [dE/dx\eff 
is the energy deposited that ends up fueling the spike 
(here nR'^yiNEexc) and m and p are close to one. They 
gave C w 0.18 for an LJ solid, which isdso appeared to 
apply to results for other pair potentialslJ. Here we use a 



picture of the ejection attained from the 2D fluid dynam- 
ics model to establish the theoretical basis for the value 
of C. That is, the internal pressure in the spike deter- 
mines a critical radius {Rg ~ Rcyi{Eexc/UoY^^) and a 
depth Rcyi, leading to the ejection of a conical volume 
of material Y sa ^RcyinRlN. This gives C = i, which is 
larger than the MD result. The difference is due in part 
to the fact that the material properties are not exactly 
those of the LJ solid and transport along crystal axes 
removes energy from the spike as discussed, however, all 
the principal features of the transport and ejection are 
the same. ThiSpCiew model resembles that of Yamamura 
and co-workersEj but disagrees with the 'so-calied' pres- 
sure pulse model used for molecular materialsES. 

Several points need further investigation. The forma- 
tion of craters at normal incidence is a topical problem 
that can be addressed by the model developed here. Fur- 
ther, the connection between the sputtering yield and 
the time used to heat the spike needs to be studied. In 
this paper, as in most of the MD simulations, we as- 
sumed it to be negligibly small. This may be correct for 
spike formation by a collision cascade, but iS|-|known to 
fail for electronic sputtering of rare-gas solidsQ. Finally, 
it must be noted that the fluid dynamics description of 
the spike is a useful complement to MD. In the fluid 
model a broad range of material properties and types 
can be readily studied, whereas complicated potentials 
are needed in MD calculations of different materials. In 
fact it is seen in Figs |[]| that the saturation leading to 
the linear regime is not simply dependent on the cohesive 
energy {E^xc ~ Uq) and the initial Rcyi, as found in the 
MD simulations using pair potentials, but also depends 
on the material parameters A and 77*. In addition, lo- 
cal equilibrium chemistry, which can play an important 
role in many of the materials of interest to us, can be 
readily included in the fluid models. However, MD has 
the advantage that non-normal incidence can be treated 
easily, the state of the ejecta (clusters vs. atoms) can 
be studied, and non-equilibrium chemistry can be intro- 
duced. Therefore, a program in which complementary 
calculations using fluid dynamics and MD simulations is 
underway. Here we have shown that a new linear sputter- 
ing regime is seen in both models and we have developed 
a simple analytic model for the yield at normal incidence. 



ACKNOWLEDGMENTS 

Part of this work was carried out during a visit per- 
formed by one of the authors (MMJ) to the School of 
Engineering and Applied Science, University of Virginia. 
Financial aid from the Astronomy and Chemistry Divi- 
sions of the National Science Foundation (U.S.A.) and 
the Consejeria de Eduacion, Cultura y Deportes del Go- 
bierno Autonomo de Canarias (Spain) are acknowledged. 



8 



* Corresponding author, e-mail: mmateo@uU.es. 

^ See for example C.T. Reimann. Mat.-fys. Medd 43, Det. 

Kong. Dan. Vid. Selsk. 351 (1993) and references therein. 
^ M. W. Thompson. Phil.Mag. 18, 377 (1968). 
^P. Sigmund. Phys.Rev. 184, 383 (1969). 
*R.E. Johnson and R. Evatt. Radiation Effects, 52, 187 

(1980) . 

^ P. Sigmund and C. Claussen. J. Appl. Phys. 52(2), 990 

(1981) . 

® H.H. Andersen. Phy.Rev.Letters 80, 5433 (1999). 
R.E. Johnson and J. Schou Mat.-fys. Medd 43, Det. Kong. 
Dan. Vid. Selsk. (1993). 

* H.M. Urbassek, H. Kafemann and R.E. Johnson. Phys. 
Rev. B 49, 786 (1994-11). 

^ E.M. Bringa and R.E. Johnson. Nucl. Instr. Methods in 
Phys. Res. B 152, 167 (1999). 

M.M. Jakas and E.M. Bringa. Phys. Rev. B 62, 824 (2000). 
M.M. Jakas. Radiation Effects and Defects in Solids, 152, 
157 (2000). 

L.D.D. Landau. Fluids Mechanics 2nd Ed., Butterworth- 
Heinemann (1995). 

Yu. B. Zel'dovich and Yu. P. Raizer. Physics of shock waves. 
Academic Press (1969). 
M.M. Jakas. In preparation. 

E.M. Bringa, M.M. Jakas and R.E. Johnson. Nucl. Instr. 



Methods in Phys. Res. B 164-165, 762 (2000). 
J.O. Hisrchfelder, C.F. Curtiss and R.B. Bird in Molecu- 
lar Theory of Gases and Liquids John Wiley & Sons, Inc. 
(1964). 

If the temperature of the spike is high, then, the speed of 
sound (c) is greater than co since c = {p/MNY^^, therefore 



according to Eq.(|2.1l[), c= (XkBT + clY^^. 



' It must be noticed that, except for a numerical factor. 



Eq.(3.2) coincides with the effective sputtering radius de- 
rived by H. Urbassek and P.Sigmund [Appl. Phys. A 33, 19 
(1984)] using a Gaussian thermal spike. As they were ob- 
tained using different models such an agreement appears to 
be a remarkable coincidence for which we have no feasible 
explanation. 

' Y. Kitazoe, N. Hiraoka and Y. Yamamura. Surface Science 

111, 381-394 (1981). 
' Z. Insepov, R. Manory, J. Matsuo and I. Yamada. Phys. 

Rev. B 61, 8744 (2000). 

H.M. Urbassek and J. Michl. Nucl. Instr. and Methods in 
Phys. Res. B22, 480 (1987). 

' D. Fenyo and R.E. Johnson. Phys. Rev. B 46, 5090 (1992- 
I) 

' E.M. Bringa, R.E.Johnson and M.M. Jakas. Phys. Rev. B 
60, 15107 (1999) 



9 



